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We investigate the nonlinear dynamics of turbulent shear flows, with and without rota¬ 
tion, in the context of a simple but physically motivated closure of the equation governing 
the evolution of the Reynolds stress tensor. We show that the model naturally accounts 
for some familiar phenomena in parallel shear flows, such as the subcritical transition to 
turbulence at a finite Reynolds number and the occurrence of a universal velocity profile 
close to a wall at large Reynolds number. For rotating shear flows we find that, depending 
on the Rayleigh discriminant of the system, the model predicts either linear instability or 
nonlinear instability or complete stability as the Reynolds number is increased to large 
values. We investigate the properties of Couette-Taylor flows for varying inner and outer 
cylinder rotation rates and identify the region of linear instability (similar to Taylor’s), 
as well as regions of finite-amplitude instability qualitatively compatible with recent ex¬ 
periments. We also discuss quantitative predictions of the model in comparison with a 
range of experimental torque measurements. Finally, we consider the relevance of this 
work to the question of the hydrodynamic stability of astrophysical accretion discs. 


1. Introduction 

Turbulent motion in the flow of an incompressible fluid has consistently eluded a 
satisfying mathematical description in over a century of investigation. The governing 
Navier-Stokes equations, although involving only a relatively simple nonlinearity, conceal 
a remarkable wealth of complex behaviour. Even in the simplest problem of the onset 
of turbulent motion in a parallel shear flow, the classical approach based on a linear 
analysis of normal modes fails spectacularly to account for the basic experimental results 
(e.g. Drazin & Reid 1981). A better insight into the transition to turbulence has been 
gained more recently through studies of the transient amplification of disturbances in 
linear theory (Butler & Farrell 1992), which, when coupled with an appropriate nonlinear 
feedback, allows perturbations to be sustained (e.g. Baggett & Trefethen 1997). The 
direct computation of nonlinear disturbances in the form of (possibly unstable) steady 
solutions or travelling waves that act as precursors of the turbulent dynamics (Nagata 
1990; Waleffe 1997) has also shed light on the process of transition. 

By its nature, fully developed turbulence demands a statistical description. Reynolds 
(1895) established the principles of a statistical theory, showing that correlations be¬ 
tween components of the fluctuating velocity field provide a stress that influences the 
bulk motion. The analogy between this turbulent transport of momentum and the vis¬ 
cous transport associated with thermal molecular motion in kinetic theory suggested the 
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concept of eddy viscosity, introduced by Boussinesq (1877). Later, Prandtl’s (1925) the¬ 
ory of the mixing length gave a predictive expression for the eddy viscosity that provided 
a remarkable quantitative agreement with experimental data such as the mean flow rate 
in turbulent pipe flow (Schlichting 1979). 

Despite the success of mixing-length theory, much effort has been expended in a search 
for more accurate representations of the Reynolds stress in turbulent flows, especially in 
engineering applications. As reviewed by Speziale (1991), the more successful approaches 
start from the exact equation governing the evolution of the Reynolds stress and apply 
a procedure of closure modelling to deal with the numerous intractable terms that arise. 
Through successive algebraic development accompanied by a large number of parameters, 
such models are able to fit an increasing range of experimental or numerical results. 
However, in our view, one disadvantage of this approach in its current state is a loss of 
physical interpretation of those terms derived from complicated algebraic constructions. 
Moreover, as noted by Speziale (1991), some of these models tend to perform poorly in 
situations for which they were not calibrated, such as rotating shear flows. 

Astrophysics provides examples of naturally occurring shear flows in which rotation 
is an essential feature. Accretion discs (e.g. Pringle 1981) are usually thin discs of gas 
in circular orbital motion around a central star or black hole, and are involved in the 
processes of star and planet formation as well as being responsible for some of the most 
luminous sources in the Universe. According to Kepler’s third law, the angular velocity of 
the gas depends on the distance from the centre as D oc and an outstanding ques¬ 

tion of potentially profound significance is whether hydrodynamic turbulence occurs in 
these situations (e.g. Balbus & Hawley 1998). While it is true that the Reynolds number 
is exceedingly large {Re > 10^^; Frank, King & Raine 2002), nevertheless the centrifugal 
instability of Rayleigh (1917) does not occur because the specific angular momentum 
increases outwards (the same criterion was derived by Solberg for homentropic com¬ 
pressible fluids; see Tassoul 1978), and no other suitable hydrodynamic instability has 
been identified. 

Historical experiments by Taylor (1923, 1936o,&) and Wendt (1933) on Couette-Taylor 
flow between differentially rotating cylinders have been adduced in the hope of a resolu¬ 
tion of this issue (Richard & Zahn 1999). The experiments suggest that turbulence can 
be sustained even in certain apparently Rayleigh-stable situations such as a Couette- 
Taylor flow in which only the outer cylinder rotates. In this case the turbulent state is 
presumably accessed through a nonlinear shear instability of the laminar state associ¬ 
ated with a subcritical bifurcation, but the wider applicability of this finding is not well 
understood. The situation is not helped by the dearth of modern experiments and the 
fact that Taylor’s findings have been challenged by Schultz-Grunow (1959). 

In this paper we investigate the nonlinear dynamics of turbulent shear flows, with 
and without rotation, in the context of a simple but physically motivated closure of 
the Reynolds-stress equation. Our approach differs from that of the conventional clo¬ 
sure models used in engineering applications. We aim to study a minimal system in 
which the modelled nonlinear terms have a clear interpretation and are as few in number 
as is compatible with the physical requirements. Indeed, in astrophysical applications 
the added complexity of other physical processes (convection, magnetic fields, radiative 
transfer, etc.) forbids anything but a minimal approach in turbulence modelling. In the 
present context of purely hydrodynamic turbulence this approach allows us to explore 
the dynamical and nonlinear behaviour in some detail without losing sight of the physical 
problem. 

The remainder of this paper is organized as follows. We formulate the model in Sec¬ 
tion 2. Next, in Section 3, we examine the local properties of the model in the contexts 
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of homogeneous shear turbulence, with and without rotation, and turbulent shear flow 
past a wall. We also compare the model with others available in the literature. Section 4 
concerns the problem of Couette-Taylor flow. The results are summarized and discussed 
in Section 5. 


2. A simple Reynolds-stress model for turbulent shear flows 

The flow of an incompressible fluid is governed by the Navier-Stokes equations 

{Of “t“ UjOj^Uf — OfpvOjjUf^ (2.1) 

^^u, = 0 , ( 2 . 2 ) 

where Ui is the velocity, p is the modified pressure (being the pressure divided by the 

uniform density p, plus the gravitational potential), v is the uniform kinematic viscosity, 
and we make use of the Cartesian tensor notation. Following a standard technique, the 
velocity and pressure may be separated into mean and fluctuating parts, e.g. 

Ui = Ui + u[, {u[) = 0, (2.3) 

where the angle brackets denote a suitable averaging procedure. We readily obtain the 
averaged Navier-Stokes equations, 

{Of + Ujdj)ui = -dip + vdjjUi - OjRij, (2.4) 

O^u, = 0, (2.5) 

where 

R^J = {uWj) ( 2 . 6 ) 

is the Reynolds-stress tensor divided by the density. 

From the fluctuating parts of the Navier-Stokes equations it is possible to obtain an 
exact equation for Rij in the form 

{Of “t“ Uj^OfflRij -\- RifzOj^Uj -\- RjkOkni vOkkRij 

= -2v{{dku'i){dkub)) - (w'Mfedfeu' -k Uju'^.dku'i) - {u{djp' + UjOtp'). (2.7) 

There is no difficulty in retaining the exact form of the linear terms on the left-hand 
side of this equation, which represent the advection of the turbulent fluctuations by the 
mean flow, their interaction with the mean velocity gradient and the viscous diffusion 
of the Reynolds stress. The difficult terms on the right-hand side cannot be written 
exactly in terms of Rij, unless further information is known about the turbulence, and 
therefore require a closure model. However, the physical effects of these terms are quite 
well understood and this insight can be used as a guide in constructing the model. 
In particular, the viscous term on the right-hand side is negative definite and causes 
a dissipation of the turbulent kinetic energy at a rate that is usually considered to be 
independent of ly in the limit of large Reynolds number. The other terms are conservative 
but allow for a redistribution of energy among the components of Rij, and it is well 
established that anisotropic turbulence has a tendency to return to isotropy (e.g. Rotta 
1951). 

Recently, one of us proposed a simple model of the stresses in astrophysical magneto¬ 
hydrodynamic turbulence (Ogilvie 2003). In the special case of hydrodynamic turbulence 
in an incompressible fluid, the model reduces to the simpler form 

{Of + Ukdk)Rij + RikOkUj + RjkOkUi = - ^R^^'^Rij - ^R^^^{Rij - ^RS^j), (2.8) 

1j Lj 
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where R = Ru is the mean-square turbulent velocity, Ci and C 2 are positive dimen¬ 
sionless constants of order unity, and L is a characteristic length-scale related to the 
geometrical constraints that limit the size of coherent structures. This model was in¬ 
tended to represent the typical astrophysical situation in which the Reynolds number is 
exceedingly large and there are no solid surfaces on which boundary layers may form. 
The modelled nonlinear terms represent two well known physical processes that are es¬ 
sential in the dynamics of turbulent shear flows. The Ci term represents the viscous 
dissipation of turbulent motion, at a rate related to the characteristic time-scale 
of the largest eddies. The C 2 term, which only redistributes energy among the compo¬ 
nents of Rij , represents the tendency of the turbulence to return to isotropy on a similar 
time-scale. This formulation gives arguably the simplest nonlinear model involving these 
two essential effects, which also guarantees that the Reynolds tensor remains positive 
definite and therefore realizable by a genuine velocity held. 

In the present paper we are interested in applying the model to laboratory shear hows 
in which the Reynolds number is not exceedingly large and in which boundary layers 
are present. We therefore enhance the model by retaining the viscous diffusion of the 
Reynolds stress and including an additional term to model viscous dissipation on the 
right-hand side: 


{dt Uk^k^Rij T Rik^k^j Rjk^kUi udkkRij 



(2.9) 


Here is a third positive dimensionless constant of order unity, and the term allows 
for the fact that, at low or moderate Reynolds numbers, when an efficient turbulent 
cascade does not form, the viscous dissipation rate is directly proportional to the viscosity. 
Hence in what follows, although we often denote as ‘turbulent’ any how for which R > 0, 
the Reynolds stresses for relatively low Reynolds numbers are more likely to represent 
the average behaviour of large-scale coherent structures (such as Taylor vortices in the 
case of the Couette-Taylor system). Wavelike behaviour, on the other hand, cannot be 
well represented in this formalism owing to the assumed locality of the dissipation. 

In the original model L was related to the thickness of an accretion disc; in a stratihed 
atmosphere it might be related to the density scale-height. It was not necessary to give 
a precise dehnition of L owing to the invariance of the original model under a rescaling 


( 2 . 10 ) 


L I—> AT, Cl I—^ AC’]^, C 2 I—^ XC2- 


In the present paper, however, we will make a definite choice for L appropriate to the 
geometrical constraints of the problem, and we will attempt to fix the values of the 
coefficients by comparison with experimental results. 

It is straightforward, as in Ogilvie (2003), to allow for a uniform rotation of the frame 
of reference with angular velocity In this case we obtain additional Coriolis terms of 
the form 


-\- Ukdk^Rij “t" Rik^k'^j T Rjk^k^i T ‘^^jkl^k^il T ‘^^ikl^k^jl “t“ * ' ' . (^■^^) 


The expression for the Reynolds-stress equation in a general orthogonal curvilinear 
coordinate system is given in Appendix A, together with an explicit expression for dkkRij 
in cylindrical polar coordinates. 
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3. Local properties of the model 

3.1. Homogeneous shear turbulence 
We first consider the idealized situation of a uniform shear flow, 


u = Sye^ 


(3.1) 


where {x, y, z) are Cartesian coordinates and the shear rate S is prescribed. We also allow 
for a uniform rotation of the frame of reference with angular velocity 


fl = Qe, 


(3.2) 


We assume that the geometrical constraints are such that L is constant. This can be 
achieved by performing a numerical simulation in a periodic box with no solid boundaries 
(Rogallo 1981; Pumir 1996), in which case L is related to the size of the box. This 
system is the only one in which the size of the coherent structures can be limited without 
imposing additional boundary conditions that would result in the loss of the large-scale 
homogeneity of the flow. The turbulence may therefore be assumed to be statistically 
homogeneous, although it is anisotropic. In this case the Reynolds stress depends only 
on time and the model reduces to a system of ordinary differential equations constituting 
an autonomous nonlinear dynamical system. We find, in detail. 


dtRccoo + 2{S -2n)R,y = 

dtRxy + 2TlRxx [S — 2Tl)Ryy = 
dtRxz + {S- 2n)Ry, = 

^tRyy “b ^^Rxy — 
dtRyz + 2TIRxz = 
dtRzz = 


(Cl -I- C2) nl/2n I £2 R3/2 _ n 

it itxx ' n T ^ r 9 t 

L 6L 

(C1+C2) i/2p Tt 

JX J^xy : 

(Cl -b C2) 1/2 „ 

ri itxz 

(C1+C2) i/2p , C2 3/2 
{Cl + C2) pl/2 p CyV 

^ n ityz ^2 

_ (Ci+C2) ^l/2^^^ ^ g^3/2 _ 


'^zz • 


(3.3) 


It can easily be seen that the two components Rxz and Ryz are decoupled from the others; 
these quantities may be expected to vanish on grounds of symmetry, as is confirmed 
below. 

The system is characterized by a Reynolds number 


V 

(3.4) 

and an inverse Rossby number 


„ -1 

(3.5) 

The Rayleigh discriminant of the rotating shear flow is 


$ = 20(2n - S). 

(3.6) 


We recall that $ < 0 is a sufficient condition for the instability of a rotating shear flow 
in the absence of viscosity. We refer informally to the situations <I> < 0, $ = 0 and $ > 0 
as ‘Rayleigh-unstable’, ‘Rayleigh-neutral’ and ‘Rayleigh-stable’ even though the criterion 
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is limited in its strict applicability. It is convenient to define a dimensionless Rayleigh 
discriminant 

(j) = Ro-^ {Ro-^ - 1) (3.7) 

and a Taylor number 

Ta = -Re^fj). (3.8) 

The dynamical system (ESJ possesses a trivial fixed point, Rij = 0, which represents 
the laminar state. Linear analysis indicates that the trivial hxed point is unstable with 
respect to infinitesimal perturbations when the Taylor number exceeds a positive critical 
value, 

Ta > Toc = (3.9) 

Although a linear instability of this kind is possible only in Rayleigh-unstable situa¬ 
tions, we demonstrate below that turbulent states can be accessed through a nonlinear 
instability of the laminar state under a wider range of conditions. 

The dynamical system may also possess non-trivial fixed points, representing states of 
statistically steady and homogeneous turbulence in which viscous dissipation is compen¬ 
sated by an extraction of energy from the shear flow. Such states may be either stable 
or unstable; even if it does not represent a statistical endpoint of the dynamics, an un¬ 
stable solution may play a transient role in the dynamics by providing an organizing 
structure in the dynamical phase space. Searching for non-trivial fixed points, we obtain 
the condition 


3L 




(Cl + C 2 ) 


rCT 


C^iy 

L2 


■ 8n(2fi - S) 


Cl ni/2 , C^iy 


(3.10) 


which may be written as a cubic equation for the dimensionless rms turbulent velocity 

u = rC^/L\S\, 



(Cl + C2)u + 




Ciu + 


a 

Re 


(3.11) 


The behaviour of it as a function of Re depends on the dimensionless Rayleigh discrim¬ 
inant (j) of the rotating shear flow. In the present model two values of (j) with special 
significance are 


C 2 


12(Ci + C2)’ 




QCi 


(3.12) 


Excluding degenerate intermediate cases, we identify four intervals of interest and illus¬ 
trate in Figure n the corresponding bifurcation diagrams in which u is plotted against 
Re. Although the qualitative features of the set of bifurcation diagrams do not depend on 
the parameters Ci , C 2 and C ^,, we make a particular selection of ‘standard’ parameters 
which is explained in Section El below. 

(«) — I <<('<<(>-■ As Re is increased, the laminar state loses stability to a branch of 
turbulent solutions at a supercritical bifurcation. 

(h) (j)- < (j) < 0. The turbulent branch bifurcates subcritically from the laminar state 
at the point of linear instability. There is an interval of Re in which stable laminar and 
turbulent solutions coexist. 

(c) 0 < (j) < (j)+. The laminar state is linearly stable for all Re and the turbulent branch 
is disconnected from it. Only the upper turbulent branch is stable, but the unstable 
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Re Re 


Figure 1. Bifurcation diagrams for homogeneous shear flow with standard model parameters 
(Cl = 0.412, C 2 ~ 0.6 and C„ = 12.48) and four different values of the dimensionless Rayleigh 
discriminant (f>. For these parameter values, cj)- — —0.049 and (j>+ — 0.243. Solid and dotted 
lines indicate stable and unstable branches. 


lower branch assists in diminishing the basin of attraction of the laminar state as Re is 
increased. 

(d) (j)+ < (j> < 00 . No turbulent solution exists and the laminar state is stable for all 
Re. 

For non-rotating shear flows, or more generally, flows with zero angular-momentum 
gradient (</> = 0), the laminar state is linearly stable for all finite values of Re. However, 
it is unstable with respect to algebraically growing disturbances at Re = 00 and a branch 
of turbulent solutions bifurcates subcritically at this point (Figure |21). 

The conclusion of this analysis is that, according to our model, statistically steady and 
homogeneous turbulence can be sustained in a rotating uniform shear flow at sufficiently 
large Reynolds number provided that the flow is Rayleigh-unstable, Rayleigh-neutral 
or else Rayleigh-stable by a sufficiently small margin. For Rayleigh-neutral or slightly 
Rayleigh-stable flows the transition to turbulence occurs through a nonlinear instability 
of the laminar state, which has a diminishing basin of attraction as Re 00 . These 
properties are in accord with the generally accepted model of transition to turbulence 
in shear flows (e.g. Grossmann 2000). Even though our model is designed principally 
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Figure 2. Bifurcation diagram for homogeneous shear flow with standard model parameters 
and vanishing Rayleigh discriminant, corresponding to a non-rotating shear flow. The branch of 
turbulent solutions bifurcates subcritically from the laminar state at Re = oo. 


to describe the statistical properties of fully developed turbulence, it appears to give a 
description of the onset of turbulence that is at least qualitatively reasonable. 

Through straightforward algebra it can be shown that the turbulent solutions share 
the following properties: 

(a) Rij is positive definite and therefore the solutions are realizable; 

(b) sigTi(—Rxy) = sign(S') and therefore the turbulent transport of momentum has the 
same sense as viscous transport; 

(c) Rxz = Ryz = 0 as expected on grounds of symmetry; 

(d) the solutions are stable with respect to arbitrary perturbations of Rxz and Ryz- 
In the limit Re —> oo equation fTTITl has at most one solution for which u tends to a 

positive limiting value. This solution exists when (j) < i.e. when the flow is Rayleigh- 
unstable, Rayleigh-neutral, or else Rayleigh-stable by a sufficiently small margin (Ogilvie 
2003). In detail, the limiting solution is 


with 


Rxx 

Rxy 

^yy 

Rzz 


3(l-Ro-^)Ci +02 


C 1 +C 2 


R 
"3 ’ 


Cl 

2LS 




/3Ro-^Ci+C2\ 

V Cl + C 2 ) 


R 
"3 ’ 


C 2 \ 

Ci + C2y 


R 
"3 ’ 


(3.13) 



6Ro"^(i?o"^ 


l)Ci 


2 c2 




4 ’) r 2 c 2 

(Cl + C2)2 


Ci(Ci+C2)2 


(3.14) 
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For a fixed Rossby number satisfying the condition (j) < (/>+, the turbulent momentum 
transport has the dependence 

- cx (3.15) 

as also occurs in Prandtl’s mixing-length theory. In our model, however, this result is 
obtained only in the limit of large Re, and the coefficient of proportionality depends on 
the Rossby number if the flow is rotating. 

The turbulent states are always anisotropic owing to the effects of shear and rotation. 
In our model, the dimensionless anisotropy tensor 

b^J = ^- (3.16) 

which describes the shape of the Reynolds tensor, depends only on the ratio CijCx and on 
the Rossby number, when the Reynolds number is sufficiently large. When CijCi is small, 
the tendency to return to isotropy is weak and the stress becomes highly anisotropic. In 
principle, the ratio CijCx could be constrained through a comparison with experimental 
data on anisotropy. We return to this point in Section rOI below. 


3.2. Turbulent shear flow past a wall 

In this section we analyse the simplest problem involving a wall-bounded turbulent shear 
flow. The solution will serve later as an asymptotic description of the turbulent boundary 
layers found in more complicated situations. 

We consider the non-rotating parallel shear flow u = Ux (y) e-x in the semi-infinite 
region y > 0 bounded by a smooth, stationary wall at y = 0 and forced by a shear 
stress Txy > 0 at j/ = oo. Even in a situation such as Couette-Taylor flow, it is usually 
permissible to neglect rotation in the turbulent boundary layers because the local shear 
rate is much larger than the rotation rate. Since the presence of the wall provides the 
only geometrical constraint on the turbulent structures, we set L = y, as is common 
in applications of mixing-length theory to wall-bounded flows. Indeed, throughout the 
remainder of this paper, we set L equal to the distance to the nearest wall. 

We seek steady solutions of the averaged equations in which the mean quantities depend 
only on y, and with the expected symmetry property Rxz = Ryz = 0. The x-component 
of the averaged Navier-Stokes equation. 


0 r'dyyUx b)yRxy , 


implies 


vdyUx — Rxy = Txy = constant. 


The non-trivial components of the Reynolds-stress equation in our model are 

op P _ (^1 + C'2) nl/2 n ^*2 3/2 n 

^-^xy^y — L ^ ^xx ' 3^^ ' ^^yy^xx jj2 ^xx') 

I ,,D ^ D 

' t'C/yyJtxy ^xyi 

I , , D ^ D 

^ vy ^ ^2 ^^yy’ 

Q ^ _ (Cl+C2) ^l/2^^^ ^ g^3/2 ^ ^dyyR^z - ^Rzz, 

subject to the no-slip boundary conditions Ux = Rij = 0 at ?/ = 0. 


R rt V - (^l+^2)ni/2n 

^yy^y^x — L ^ ^xy 

Q ^ _ (Cl+C2) ^l/2 ^ ^^3/2 


(3.17) 

(3.18) 


(3.19) 
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We rewrite the equations in a dimensionless form by means of the standard transfor¬ 
mations 

ux{y) = v{r])y^, R^j{y) = nj{r])Txy, tj = , ( 3 . 20 ) 

Observing the property = ryy, which implies r = r^x + ‘^fyy, we obtain the reduced 
problem 


- 


TvvV = -- 


0 = -- 


V - Vxy 

= 

1 , 



*1/2^ 

' ' XX 

L 

^2 ^3/2 

+ r" - 

' XX 

c. 


3?y 

7y2 

•^1/2 
' ' xy 



+ r" - 

' ' xy 

c. 



9 

7?^ 



^2 ^3/2 

1 

c. 

' 'yy 


3?y 

^vv 

9 ' 

jjZ 


(3.21) 


(3.22) 


with boundary conditions u( 0 ) = 'ey ( 0 ) = 0 . 

The desired behaviour at 77 = 00 can be deduced by analysing the far-held limit ry ^ 1, 
for which we obtain the asymptotic form 


v'^T] ^ 


-V 2 II ^ 


with 


Gj - 

I'ijO + Cj 

iV ^ 4 

/ Cl 3/2 

ro=( 

' 6 

2 > 

C 1 C 2 

_(3Ci+C 2)^ 
3(Ci-kC2)'’°’ 

^xyO — 

-1, 


1/2 


(O1+C2), 


C 2 


ryyO — 


r^o- 


(3.23) 


(3.24) 


(3.25) 


3(01 + C 2 ) 

The reduced problem is universal, involving no parameters other than the model pa¬ 
rameters Cl, C 2 and Cu. To solve it numerically, we use a hnite-difference Newton- 
Raphson relaxation method on a stretched mesh, starting from an initial guess with 
rij = rijo and v = rj. The outer boundary condition rb(77out) = 0, where ryout ^ 1, 
imposes the desired far-held behaviour with adequate hdelity. We choose ryout = 10^. 
The desired universal boundary-layer solution is shown in Figure 13 The stability of this 
solution has been conhrmed using a time-dependent numerical method. 

In the Prandtl-von Karman analysis of turbulent boundary layers (e.g. Schlichting 
1979) the velocity prohle for ry ^ 1 is given as 


V Alnr] + B 


(3.26) 


where A and B are dimensionless empirical constants, and it is customary to refer to k = 
\/A as the von Karman constant. This ‘universal velocity prohle’, ‘log law’ or ‘law of the 
wall’ is generally in excellent agreement with experimental data. The values traditionally 
assigned on the basis of Nikuradse’s experiments in the 1930s are k = 0.41 and B = 5.2 
(for a smooth wall). However, recent experiments at higher Reynolds numbers show a 
much superior ht with k = 0.436 and B — 6.15 (Zagarola & Smits 1998). 

In our model the integrated velocity prohle for 7 y ^ 1 is 

v{'q) = Vo + v[\TLy- —+ 0{iri~‘^), (3.27) 

which is asymptotically equivalent to the Prandtl-von Karman velocity prohle; thus we 
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Figure 3. Universal boundary-layer solution for turbulent shear flow past a wall with standard 
model parameters. The dimensionless Reynolds stress components r^x, Tyy = r^^ and —r^y are 
shown. 

identify the von Karman constant as 

1 2 

The von Karman constant depends only on Ci and C 2 and can therefore be fitted inde¬ 
pendently of Cu- Figure^shows the relation between Ci and C 2 for constant k. 

The additive constant B = vq cannot be deduced from the asymptotic analysis but 
must instead be determined from a numerical solution of the problem including the 
viscous sublayer close to the wall. Numerical integration of the boundary layer equations 
shows that vq depends primarily on Cjy, as shown in FigureEl For simplicity, in that figure 
we have reduced the parameter space by imposing the constraint n = 0.436 (Zagarola & 
Smits 1998). 


C1C2 


Q{Ci + C2fl 


(3.28) 


3.3. Comparison with other models 

In comparison with other Reynolds-stress models available in the literature, our model 
appears quite simplistic (e.g. Speziale 1991). The viscous dissipation rate is given by 


2L 2L2 


R 


(3.29) 


and we do not attempt to model a separate time-dependent equation for this quantity. 
One reason for this is that the length-scale L is imposed by the geometrical constraints 
in our problem, and is not free to expand as occurs when turbulence is generated in 
a localized region within a larger system. Instead, we allow for the effects of a finite 
Reynolds number by specifying a dissipation time-scale that is related to the character¬ 
istic time-scale of the largest eddies in high-Reynolds-number turbulence, and to 

the viscous time-scale /v at lower Reynolds numbers. 

In freely decaying turbulence with no mean shear or rotation, the return to isotropy is 
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Figure 4. Contours of constant k with (from bottom to top) k = 0.5, 0.46, 0.436 (solid line), 

0.42, 0.38, 0.34 and 0.3. 
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Figure 5. Contours of constant vo as a function of CijCi (for a fixed k = 0.436) and \/C7. 
The solid line corresponds to the experimentally determined value of uo = 6.15. The uniform 
contour spacings suggest that vo depends primarily and approximately linearly on \/C7. 


described in our model by the equation 


dbij R dbij 


( 2C2 ^ , 

VCi + auL-ii?-i/2 ) 


dr e dt 


(3.30) 
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where 



(3.31) 


is the anisotropy tensor and r is a dimensionless time variable (e.g. Speziale 1991). In the 
limit of large Reynolds number we therefore obtain a linear return to isotropy, equivalent 
to the one introduced first by Rotta (1951). 

Furthermore, we do not attempt to model the ‘rapid pressure-strain correlation’, for 
which elaborate algebraic models have been proposed (e.g. Sjogren & Johansson 2000). 
Through this simplihcation we may lose some accuracy. However, it is the treatment of 
this term that has given rise to models that deal poorly with rotating shear flows. For 
example, the widely adopted model of Launder, Reece & Rodi (1975) is not consistent 
with Rayleigh’s criterion in the sense that it does not permit turbulence to be sustained 
at high Reynolds number in certain situations where Rayleigh’s stability criterion is 
not satisfied (Speziale 1991), perhaps because a realizability condition of some kind is 
implicitly violated. It is to be hoped that a way can be found to model the rapid pressure- 
strain correlation in future with due regard to Rayleigh’s criterion. 


3.4. Summary of the constraints on the parameters Ci, C 2 and C^, 


We have purposely adopted a simple closure model for the Reynolds-stress equation 
so that we can focus on the nonlinear dynamical properties of the system rather than 
engaging in a lengthy exercise of parameter fitting. Nevertheless, because the model 
naturally predicts a logarithmic velocity profile close to a wall, it makes sense to apply 
the two accurate constraints k = 0.436 and vq = 6.15 provided by the very high-quality 
experimental data on wall-bounded turbulent shear flows in the Superpipe experiment 
(Zagarola & Smits 1998). The first constraint provides a relation between Ci and C 2 only, 
whereas the second (which applies only for a smooth wall) yields C^, provided that Ci and 
C 2 are known. A third constraint, which would be required to fix all three parameters 
of our model, might in principle be provided by experimental data on the anisotropy 
of the Reynolds stress in homogeneous shear turbulence, or on the return to isotropy of 
homogeneous turbulence. In fact, the limitations of the three-parameter model mean that 
no choice of the parameters can accurately match all experimental results. For example, 
Choi & Lumley (2001) find that the return to isotropy of homogeneous turbulence is 
more complicated than is assumed in any available closure model. In Section 14.21 below 
we compare the predictions of our model with data from Couette-Taylor experiments, 
and tentatively deduce an approximate value of C 2 — 0.6. Hence in what follows (unless 
otherwise mentioned) we shall take as standard parameters Ci — 0.412, C 2 = 0.6 and 
Cy = 12.48, which yield k = 0.436 and Uq = 6.15 as required. The predicted boundary- 
layer velocity profile with this choice of parameters is compared with the experimental 
measurements of Zagarola & Smits (1998) in Figure |S| 


4. Couette-Taylor flow 

Couette-Taylor flow between differentially rotating coaxial cylinders is a seemingly 
simple dynamical system that has been found to exhibit a rich variety of nonlinear 
behaviour. Much of this interesting dynamics occurs close to the onset of Rayleigh’s cen¬ 
trifugal instability, albeit in a confined setting and in the presence of viscosity. Our main 
interest here is in the existence and properties of turbulent states in Couette-Taylor flow 
at large Reynolds numbers, rather than the behaviour close to the onset of instability. 
This aspect of the problem has received rather little attention from experimentalists or 
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Figure 6. Solution of the dimensionless boundary-layer problem for turbulent shear flow past 
a wall with standard model parameters. The dotted line corresponds to the asymptotic profile 
V = (0.436)“^ ln77 -b 0.615 (Zagarola & Smits 1998) which is their suggested best fit for their 
experimental results. The small black dots are the experimental data points from Zagarola & 
Smits, corrected for systematic errors according to the prescriptions of McKeon et al. (2003). 
The open triangles are the experimental data points from Reichardt (1940). 

turbulence modellers. Recently, however, arguments based on the Couette-Taylor sys¬ 
tem have been made in connection with important questions relating to turbulence in 
astrophysical flows involving differential rotation (e.g. Richard & Zahn 1999). 

Since our Reynolds-stress model is based on a covariant formulation, it naturally in¬ 
cludes the effect of rotation as well as shear. In this section we apply the model to the 
Couette-Taylor system, compare its predictions with the available experimental results 
and examine the wider consequences of these findings. 

4.1. Predictions of the model 
4.1.1. Governing equations and numerical solution 

We consider the shear flow between two infinite coaxial cylinders located at radii ri 
and To, rotating with angular velocities and Oq. Adopting a cylindrical coordinate 
system (r, 0, z), we seek steady solutions of the averaged equations in which the mean 
quantities depend only on r and the mean flow is azimuthal only: u = rQ{r)e^. In 
this case the Reynolds-stress equation in our model reduces to and straightforward 
algebra yields Rrz = R 4 .Z = 0. We choose the scale-length L to be the distance to the 
nearest wall, namely L = min(r — ri^ro — r). 

Experiments on Couette-Taylor flow, using a cylindrical container of finite height 
h, are perturbed by end-effects, in particular the Ekman circulation, when the aspect 
ratio h/{ro — ri) is not very large. We do not attempt to model end-effects, but note 
their potential significance when comparing our findings with experimental results (see 
Section lOHl . 

The angular velocity profile 11 (r) between the cylinders is obtained self-consistently by 
solving also the angular momentum conservation equation (the azimuthal component of 
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the averaged Navier-Stokes equation), 

{r"^Rr4, + vr'^S) = 0 (4.1) 

where S = —rdU/dr is the shear rate in cylindrical geometry. Note that equation 14.111 
can be integrated to introduce the torque T between the cylinders, 

r'^ + vr'^ S =. (4.2) 

zirnp 

We obtain a ninth-order system of ODEs with one eigenvalue (T) which requires ten 
boundary conditions: Rrr = Rnp = Rci>^ = Rzz = 0 on each boundary, as well as 
r2(ri) = fli and f2(ro) = flo- We solve this two-point boundary-value problem numerically 
using a Newton-Raphson relaxation method, using the solution for Re —> oo as an 
initial guess (or, whenever applicable, the results of a previous calculation for similar 
parameters). 

In addition to possible turbulent solutions, there is of course the well-known laminar 
Couette-Taylor flow fl(r) = a + , where a = {TIotI — flirf){rg — rf)~^ and (3 = 

— no)rfr^{rg — rf)~^. The local dimensionless Rayleigh discriminant of the laminar 
solution is (j){r) = {a/ -|- (a//l)r^. Unlike plane Couette flow, the laminar Couette- 
Taylor flow is linearly unstable in certain regions of the parameter space. Therefore 
turbulent states may be accessed through either linear or nonlinear instabilities. 


4.1.2. Asymptotic analysis 

The results presented in Section E3 suggest that the system of equations (Tra and 
could also be solved approximately by asymptotic matching, between the universal 
boundary-layer solution near the wall and a high-Reynolds number limiting solution 
{Re oo) in the main body of the fluid. When Re oo, the turbulent solution is 

(4.3) 

whenever this is positive, and then 


C2-6Ro-flRo-^ - l)Ci 
Ci(Ci+C2)2 


Rrd> — 


^ 3/2 _ ^1 


2LS 


^2^/2 

'C 2 - 6C'iRo-i(i?o-i - 1)' 

Uj 

[ Ci(Ci+C2)2 J 


1 3/2 




(4.4) 


by direct analogy with results and (liTTKll of the local analysis. Unfortunately, 

there exist no analytical solutions to the angular momentum equation with this Reynolds 
stress prescription unless Ro ^ 1. The Rossby number of the laminar flow is large for 
typical narrow-gap setups; one might therefore hope to use this asymptotic limit for 
the study of the turbulent regime also. However, numerical solution of equations jSD) 
and (lOl reveals that turbulence effectively reduces the shear outside the boundary 
layers and prevents the use of the Ro 3> 1 asymptotic limit unless the gap is extremely 
small (typically, less than a few percent of the average radius). For completeness, we 
nevertheless provide such an asymptotic analysis in Appendix B. 


4.1.3. Stability diagram and structure of solutions across parameter space 

In what follows, we call the Couette-Taylor flow ‘unstable’ whenever there exist solu¬ 
tions to the equations of the model with R > 0. By doing so, we implicitly assume that 
the background noise level is sufficiently high to excite finite-amplitude instabilities if 
they exist. Figured shows a stability diagram for corotating cylinders in the {ReojRefl 
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Figure 7. Stability boundaries for Couette-Taylor flow in a flxed geometrical setup 
{ ri/vo = 0.7), predicted with standard model parameters. The black symbols are the predictions 
of our model, and delimit the regions of turbulent solutions in the top-left and bottom-right cor¬ 
ners. The open symbols are the data from Richard (2001) for the same geometrical setup. The 
solid line is the stability limit according to Rayleigh’s criterion = floUo), the dashed line 

marks the Keplerian ratio = Q,oro^^) and the dotted line marks solid-body rotation 

(fli = flo). 


plane (where Rco = droflo/u and Rci = driPlijv^ with d = Tq — Vi) for a given geometri¬ 
cal setup {rijro = 0.7). As predicted by the local analysis, Rayleigh-stable flows can be 
subject to finite-amplitude instabilities provided they are notionally Rayleigh-stable by 
a small margin only, thus displacing the stability boundary in the top-left corner of the 
diagram by a small amount. Strong enough shear can overcome the stabilizing angular- 
momentum gradient in the case where Reo 2> Ret and finite-amplitude instabilities are 
also found in that region (bottom-right corner of the diagram). The predicted onset of 
instability in the case of relatively low Reynolds numbers and also of counter-rotating 
cylinders is discussed further in Section l4.2.dl Comparison with experimental data from 
Richard et al. (2001) shows significant discrepancies with the predictions of the model, 
though this could be expected because the aspect ratio of their experimental setup is not 
large (see the discussion in Section B.2.1ll . 

The structure of the solutions in various regions of parameter space, as shown in 
Figure IHl reflects the physics of turbulent flow. In the case where the outer cylinder is at 
rest, the local dimensionless Rayleigh discriminant for the laminar solution is 

(j)[r) = (r/vof - {r/rof (4.5) 

and is therefore always negative (Rayleigh-unstable). Viscosity stabilizes the flow for 
sufficiently low Reynolds number Rci but the transition to turbulence occurs directly, 
through a linear instability (the critical value for the transition depends on the geomet¬ 
rical setup). Turbulent stresses are largest in the bulk of the fluid, near the mid-point 
Tm = (ro + n)l‘^- 
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Figure 8. Structure of solutions for a fixed geometrical setup with ri/ro = 0.7, for three 
regions of parameter space in the case of corotating cylinders. The left-hand column shows 
the local dimensionless Rayleigh discriminant <?i(r); in each figure, the horizontal dotted line 
marks the position of the critical value 4>+ = Cij^Cx ~ 0.243. The right-hand column shows the 
mean-square turbulent velocity, normalized to its maximum value for clarity. The top two panels 
(a) correspond to the case where the outer cylinder is at rest, with the laminar solution marked as 
a dotted line, then the outer cylinder rotation rate is steadily increased with Reo = 10^ (dashed 
line), Rco = 10® (short-long dashed line), Rco = 10^ (long-dashed line) and Rco = 10® (solid 
line). The line-style coding is the same for both plots. Panels (b) correspond to the case when 
the inner cylinder is at rest, showing the laminar solution (dotted line), then turbulent solutions 
for Reo = 10® (dashed line), Reo = 10® (short-long-dashed line), Reo ~ 10^ (long-dashed line) 
and Reo = 10® (solid line). Panels (c) correspond to the onset of instability near the Rayleigh 
stability limit. The outer cylinder rotation is fixed {Reo = 10®) and the inner cylinder rotation 
is varied. The solution is laminar (dotted line) for Ret = 1.22 x 10®, then becomes successively 
more turbulent for Rd = 1.24 x 10® (dashed line), Rei = 1.26 x 10® (short-long-dashed line) 
and Rei = 1.28 x 10® (long-dashed line). The last curve lies further away from onset, and into 
the Rayleigh-unstable zone with Rei — 2 x 10® (solid line). 
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In the case where the inner cylinder is at rest, 

(j){r) = (r/nf - {rlnf (4.6) 

is always positive (Rayleigh-stable). Linear instability is therefore not expected, but 
turbulent states may still be accessed through a nonlinear instability. Sufficiently close 
to the inner cylinder, in the region where (j) < (/)+, the local analysis would suggest 
that the laminar solution is unstable to finite-amplitude perturbations provided the local 
Reynolds number is large enough. Depending on the gap width, two situations may arise: 
either 0 < (jj < (j)^ for all r, or there exists a transition within the fluid between a locally 
stable region and a region of finite-amplitude instability. An equivalent way of looking 
at the problem is to note that the stabilizing effect of rotation on the development of 
turbulence is weaker near the inner boundary when the inner cylinder is at rest, and so 
we expect turbulence to develop first near the inner cylinder (as can indeed be seen in 

Figure IHJ. 

Finally, we explore the behaviour of solutions near the onset of nonlinear instability at 
very high Reynolds number in the region close to the Rayleigh line (flirf = floCo). The 
transition to turbulence occurs when (j) for the laminar solution near the inner cylinder 
drops below This happens in the Rayleigh-stable domain. As Rci is increased, (j) on 
the outer cylinder drops below 0 which marks the transition to the Rayleigh-unstable 
domain. 

The structure of the solutions in all three cases seems to suggest that turbulence is 
extremely efficient in transporting the applied torque: the marginal solution (j) = (j)+ is 
favoured near onset and also in the Rayleigh-stable case far from onset. In the latter 
case the solution deviates from the marginal stability solution only in the thin viscous 
boundary layers. This behaviour is typically also observed in convection. In the Rayleigh- 
unstable case on the other hand, a solution with (j) = (j)+ >0 could not possibly satisfy 
the applied boundary conditions, and the flow appears to compromise by choosing an 
intermediate solution with (/) < 0 in a significant part of the domain. 

4.2. Gomparison with experimental data 
4.2.1. Discussion of available data 

Since Taylor’s (1923) pioneering work on the stability of fluid flows between two coax¬ 
ial rotating cylinders, a wealth of experimental data has been collected on the dynamical 
properties of such flows for various aspect parameters (vi, Tq, h) and for a very large region 
of the (Rco, Rei) parameter space. In particular, Wendt (1933) and Taylor (1936a, &) pre¬ 
sented the most extensive collection of torques and velocity measurements for turbulent 
Couette-Taylor flow far from onset, whereas Andereck, Liu & Swinney (1986) reviewed 
the successive flow-pattern transitions near onset. 

Amongst other notable results, Wendt (1933) studied the contaminating Ekman flows 
arising from end-effects, which can drive significant deviations from the laminar Couette- 
Taylor flow. He proposed an ingenious system including a free top surface and a differen¬ 
tially rotating split bottom plate to reduce end-effects; this setup is indeed able to reduce 
meridional flows but not completely. Comparisons of torque measurements between var¬ 
ious bottom boundary conditions revealed that end-effects are especially important for 
aspect ratios h/d smaller than 40. Wendt performed experiments with only the outer 
cylinder rotating, and showed that torque measurements made with bottom plates coro¬ 
tating with the outer cylinder were roughly 10% larger than in the case where the bottom 
plates are stationary for an aspect ratio of 50, 100% larger for an aspect ratio 23 and 
400% larger for an aspect ratio of 11. Naturally, any theory that assumes axial trans¬ 
lational symmetry for the system (as does the model investigated here) can only be 
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compared with experiments that have little contamination from end-effects, and we may 
use Wendt’s findings as a guideline for distinguishing between adequate and inadequate 
sets of experiments. 

Experiments at very high Reynolds number, and with a wide gap {rijro = 0.724), have 
been performed in the case where only the inner cylinder is rotating by Lathrop, Fineberg 
& Swinney (1992) and more recently by Lewis & Swinney (1999). Richard et al. (2001) 
have performed such experiments (rijro = 0.7) with both inner and outer cylinders 
rotating, using split-bottom boundary conditions. Wide-gap setups are well suited to 
verify the adequacy of our theory in describing the effects of rotation on turbulent shear 
flows (as opposed to narrow-gap setups, which are principally dominated by the shear 
instability). However, all these experiments have too small (< 25) an aspect ratio to be 
free of end-effects, so we cannot use them reliably for comparison with our work. 

In what follows fSection I4.2.2II . we first compare the predictions of our model with the 
torque measurements from Wendt (1933) in order to obtain a third constraint on our 
parameters C\^Ci and Cv, as anticipated in Section E3I We then compare the model to 
Taylor’s (1936a) data. In Section [4.2.31 we look at the predictions of the model for the 
onset of linear instability in Couette-Taylor flows and compare it with Taylor’s (1923) 
experimental work. Finally, we discuss the stability of Keplerian shear flows in the light 
of our model, and compare our predictions with those of Richard & Zahn (1999) in 
Section 


4.2.2. Torque measurements 


Wendt’s (1933) torque measurements present the most extensive results for an experi¬ 
mental setup with large aspect ratio. By using his narrow-gap data, which are relatively 
free from contaminating end-effects, we attempt to constrain our basic parameters fur¬ 
ther. Wendt presents the results of twelve sets of experiments for the following setup: 
h = 50 cm, To = 14.70 cm and = 13.75 cm. For each set of measurements, the ratio of 
the rotation rates of the inner and outer cylinders is fixed, and the Reynolds numberf is 
defined as 


Re = 


^i\r 


(4.7) 


Wendt plots in his Figure 9c the ratio of the turbulent to the laminar torque. We have 
performed ten suites of numerical calculations where C 2 is varied between 0.1 and 1 
by increments of 0.1, and for each set have calculated the typical error between the 
predictions of our model and Wendt’s experimental data points with the formula 


EiC2) = E 


T 

T 


Tlam 

mod 

exp 


(4.8) 


where the sum spans all data points in all twelve sets of experiments. The results are 
shown in FigureEland suggest that the best fit is obtained with parameters Ci = 0.412, 
C 2 = 0.6 and C^, = 12.48, which we have adopted as standard throughout this paper. The 
corresponding fit to Wendt’s experimental data is shown in Figure 11(11 We note that the 
fit is quite good though we are not able to fit all the curves equally well. In particular, the 
experimental results for a stationary inner cylinder (black circles on the right-hand panel) 
seem to deviate significantly from the predictions of our model for all possible values of 
C 2 - Leaving this particular set of experiments out of the least-square fitting procedure 
seems to reduce the optimal value of C 2 , but not significantly (see the square symbols in 
FigurelHl which have a minimum near C 2 = 0.55). We emphasize that the constraint on 

t Wendt actually uses a quantity related to the Reynolds number, (60/27r)|f2o — flil/n. 
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Figure 9. Study of the model parameters that best fit Wendt’s (1933) experimental data. For 
each value of C 2 , C\ is chosen such that k = 0.436 and Cv is chosen such that uo = 6.15. The 
corresponding values of Ci and Cv are shown as solid lines, with the relevant scales on the left 
and right of the plot respectively. Using these parameters, the error E (as given by equation 
EiU is calculated and shown as triangular symbols. The minimum error occurs for C 2 = 0.6. We 
also plot the error E calculated when leaving out the experimental results for stationary inner 
cylinder as square symbols. 


our parameters obtained by fitting Wendt’s data is much less satisfactory than the two 
tight constraints provided by the universal velocity profile of turbulent boundary layers. 
Therefore the estimate C 2 = 0.6 is to be regarded as tentative and approximate only. 

We then tested the predictions of the model against Taylor’s (1936a) torque mea¬ 
surements. Taylor’s data consists of eight sets of experiments for varying gap width; he 
compares the torques measured for similar Reynolds numbersf (as defined by equation 
when only the inner cylinder is rotating, and when only the outer cylinder is rotat¬ 
ing. The results are presented in Figure im the predictions of the model show excellent 
agreement with the experimental data, even near onset, in the case where only the in¬ 
ner cylinder is rotating. The agreement is also very good (except for the two widest gap 
widths) in the case where only the outer cylinder is rotating, except near onset. However, 
Taylor reports that the onset of instability in the case where only the outer cylinder is 
rotating undergoes hysteresis (where the turbulent solution can only be accessed through 
finite-amplitude perturbations); since our model assumes that the turbulent solution is 
chosen whenever it exists, we represent only the turbulent branch of the hysteresis loop, 
whereas it appears that Taylor’s data follows the laminar branch up to a critical Reynolds 
number that may depend on the amount of noise present in his apparatus (see Schultz- 
Grunow 1959 for an assessment of the critical Reynolds number for the persistence of 
laminar flow in a noise-free Couette-Taylor system). 


t Taylor actually uses a quantity related to the Reynolds number, H/27ru. 
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Figure 10. Comparison between the predictions of the model for Ci — 0.412, C 2 ~ 0.6 and 
Cl, = 12.48 and Wendt’s (1933) data extracted from his article, Fig. 9c. The left panel corre¬ 
sponds to his measurements for counter-rotating cylinders, and the right panel corresponds to 
corotating cylinders. The corresponding line-styles and symbols are shown in each diagram. 


4.2.3. Onset of instability in Couette-Taylor flows 

Although our model was initially designed for fully developed turbulent flows at very 
high Reynolds number, and in fact in the context of astrophysical magnetohydrodynam¬ 
ics (Ogilvie 2003), we now show that the addition of the viscous correction terms (see 
Section|21) is also able to reproduce (qualitatively, and to some extent also quantitatively) 
the onset of instability in the Couette-Taylor system. 

The classic work of Taylor (1923) combined experimental and theoretical studies of 
the onset of linear instability in the laminar Couette-Taylor flow for both corotating and 
counter-rotating cylinders. Impressive agreement was found between the appearance of 
axisymmetric Taylor vortices in the experiments and the occurrence of an axisymmetric 
linear instability in the theoretical calculation. We have investigated the linear stability 
of the laminar flow within the context of our model. To do this, we linearize the Reynolds- 
stress equation about the laminar solution and seek solutions of the form Rij = Rij (r) e®*. 
The linearized system of ordinary differential equations admits a set of discrete modes, 
with the growth rate s appearing as an eigenvalue. We solve this system numerically and 
identify the stability boundary as the position in the parameter space where the largest 
eigenvalue passes through zero. The results depend on C,,, but not on Ci or C 2 as these 
two parameters appear only in nonlinear terms. 

The linear stability boundary predicted by our model is shown in Figure IT^ in com¬ 
parison with Taylor’s experimental results. For C^, = 12.48 the agreement is quite good 
(an even better fit can be obtained for C^, = 11). We therefore again find that, although 
our model is designed principally to describe fully developed turbulence, it also performs 
quite well in describing the onset of instability. Of course, Taylor vortices themselves are 
not a turbulent flow, but our model does not make a clear distinction between coherent 
and turbulent flows at relatively low Reynolds numbers. The near coincidence between 
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Figure 11. Comparison between Taylor’s original data (1936 a) and the predictions of our model 
with standard parameters. For each panel, the outer cylinder has radius To = 4.05 cm and the 
inner cylinder radius Vi is given in cm. The upper branch corresponds to the torques measured 
on the outer cylinder when the inner cylinder is rotating, and the lower branch corresponds to 
torques measured on the inner cylinder when the outer cylinder is rotating. The angular rotation 
rate 11 is related to as 12 = 2'kN. The dotted line shows the laminar solution whereas the 
solid line is the prediction of the model. The quantity T/N^ is measured in CCS units. 


our linear stability results and Taylor’s is not trivial because, unlike him, we do not 
represent or solve for the optimal axial wavenumber of the linear disturbance. 

Also shown in Figure lT^ is the nonlinear stability boundary which delimits the region of 
parameter space in which our model predicts turbulent solutions to exist. The discrepancy 
with Taylor’s results illustrates the fact that a finite-amplitude instability, apparently 
not detected in Taylor’s (1923) experiments, may occur in the case of counter-rotating 
cylinders. This idea is supported by the results of Coles (1965), who reports on the 
existence of a well-defined hysteresis zone delimited by a boundary qualitatively similar 
to our nonlinear stability boundary. 
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Figure 12. Comparison between the predictions of our model for linear and nonlinear stability 
of Couette-Taylor flow between two cylinders with ro = 4.035 cm and rt as shown in the 
plots, and Taylor’s (1923) experimental data for the onset of instability. The black symbols 
correspond to Taylor’s data, the open symbols corresponds to the onset of nonlinear instability 
as calculated with our full model for standard parameters. The solid line corresponds to the 
onset of instability calculated through a linear stability analysis of our turbulent model for 
Cv = 12.48 and the dotted lines show the same thing for Cv = 11. 


This boundary turns over and crosses the Rci = 0 axis for finite Rco- The unstable 
domain thus delimited for Rei < 0, Rco < 0 is the point-symmetric domain to the one 
identified as a region of finite-amplitude shear instability in the quadrant Rci > 0, Rco > 
0 (see Section EH and Figure [TJ. 


4.2.4. Keplerian shear flows 

An important unsolved problem in astrophysics concerns the hydrodynamic stability of 
accretion discs in which gas flows in circular Keplerian orbits with Q oc Although 

magnetohydrodynamic instabilities are known to be effective in generating turbulent mo¬ 
tion and angular momentum transport in discs that are sufficiently ionized (e.g. Balbus 
& Hawley 1998), such mechanisms probably fail to operate in some important circum¬ 
stances such as in very weakly ionized regions of protoplanetary discs. 

Recently, Richard & Zahn (1999) suggested, not unreasonably, that the stability of 
Keplerian flows might be deduced from the results of wide-gap Couette-Taylor exper¬ 
iments. They extract from Taylor’s (1936a) and Wendt’s (1933) data that the critical 
Reynolds number for instability, in the case where the inner cylinder is at rest, varies 
roughly as for wide gaps. From this result they deduce that there must exist a 

local Reynolds number for rotating shear flows 


Rgrz = — 

V 


an 

dr 


(4.9) 


with a critical value Rec,RZ — 6.3 x 10® for instability (see Figure El- If such an ab¬ 
straction of the Couette-Taylor experiments to Keplerian flows is indeed justified, this 
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Figure 13. Critical Reynolds number Re^ = min(Reo) for the onset of nonlinear instability in 
the case where the inner cylinder is at rest, as a function of gap width {vo = 4.05 cm for all 
points, and ri is varied between 3.2 cm and 3.94 cm). The error bars reproduce Taylor’s own 
interpretation of his data (1936a) with the size of the error bar corresponding to the extent 
of the hysteresis loop. The dashed line is the fit proposed by Richard & Zahn (1999), with 
Rcc = 6.3 X 10®(d/rm)^; the solid line is our own results, and the inclined dotted line is a 
fit for the wide-gap limit with Rcc = 2 x ■ The two horizontal lines are critical 

Reynolds numbers for plane Couette flow: Rcc = 1300 is derived from Dauchot & Daviaud’s 
(1995) experiments, and Rcc = 412 is the prediction of our model. 


criterion would suggest that Keplerian discs (which typically have Reynolds numbers 
many orders of magnitude larger than this critical value) are indeed likely to be tur¬ 
bulent. Numerical solutions of our model near the onset of nonlinear instability in the 
case where the inner cylinder is at rest also reveal that the critical Reynolds number for 
instability varies as for wide gaps, although the proportionality constant is lower 

(see Figure ICTl. The same argument proposed by Richard & Zahn (1999) applied to our 
results would therefore also suggest that Keplerian shear flows should be unstable. How¬ 
ever, a local analysis of our model predicts instability for Keplerian shear flows in the 
limit of large Reynolds numbers only when C 2 /C 1 >8/3 (Ogilvie 2003), which is not the 
case for the standard model parameters chosen in this numerical experiment. Hence, it 
is not clear that a generalization between Couette-Taylor results with a stationary inner 
cylinder and Keplerian flows can be made. 

More generally, it is not clear that stability results for wall-bounded flows can be ap¬ 
plied to unbounded flows. The instability may be triggered precisely by the presence of 
the boundaries (both the side walls, through the non-local effect of a redistribution of the 
shear profile between the cylinders, and the bottom boundary, through contaminating 
Ekman flows). For instance, in Figure^lwe explore the stability of Couette-Taylor sys- 
terns with inner and outer cylinders in Keplerian ratios (fljrd = VtoVo )■ As mentioned 
earlier, if walls were absent and the shear was everywhere Keplerian, the local analysis 
of our model would only predict instability if C' 2 /C'i > 8/3. However, we find that for 
large enough Reynolds numbers, instability can be found for ratios of C 2 IC 1 < 8/3 in a 
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Figure 14. Predictions of the model for the onset of instability in a Couette-Taylor system 
identical to that of Richard (2001) with ro = 5 cm, rt = 3.5 cm. For each enter cylinder rotation 
rate (represented by its Reynolds number Rco), the curves represent the ratio of the critical 
Reynolds number Rsi^crit for instability to the Reynolds number corresponding to a Keplerian 
rotation rate for the inner cylinder Rei^kep = Reo{rolrifl^^. 


wall-bounded experiment. This behaviour is possible because the dimensionless Rayleigh 
discriminant (j) of the laminar solution is not uniform. 


5. Discussion 

In this paper we have investigated the nonlinear dynamics of turbulent shear flows, 
with and without rotation, in the context of a simple but physically motivated closure of 
the equation governing the evolution of the Reynolds stress tensor. In order to permit a 
detailed exploration of the nonlinear behaviour and to emphasize the physical interpre¬ 
tation of the dynamics, the approach we have taken differs from that of the conventional 
closure models used in engineering applications. We have not developed a model of great 
algebraic sophistication and attempted to fit the large number of parameters therein 
by applying a restricted class of constraints. Instead, we have adopted a minimal clo¬ 
sure of the Reynolds-stress equation in which the modelled nonlinear terms have a clear 
interpretation and are as few in number as is compatible with the physical requirements. 

Our model, equation lES, retains the exact form of the linear terms representing 
the advection of the turbulent fluctuations by the mean flow, their interaction with the 
mean velocity gradient and the viscous diffusion of the Reynolds stress, while using a min¬ 
imal set of algebraic terms with three dimensionless parameters to represent dissipation 
through a turbulent cascade (with parameter Ci) and through direct viscous damping 
(parameter C^), as well as the tendency to return to isotropy (parameter 02 ). 

In a local analysis of homogeneous shear turbulence with or without rotation (Sec- 
tion ld.lll . our closure model reduces to an autonomous nonlinear dynamical system whose 
fixed points, either stable or unstable, represent the laminar state and any statistically 
steady turbulent states. We find that the behaviour of the system depends on the Rayleigh 




26 


P. Garaud & G. I. Ogilvie 

discriminant (defined by equation Id.611 i of the rotating shear flow. The model predicts 
that Rayleigh-unstable flows become turbulent at sufficiently large Reynolds number 
through a linear instability associated with a supercritical (or, rarely, subcritical) bifur¬ 
cation. Flows that are Rayleigh-stable by a sufficiently large margin are predicted not 
to support sustained turbulence however large the Reynolds number. This behaviour is 
naturally consistent with Rayleigh’s stability criterion. 

Non-rotating (Rayleigh-neutral) shear flows and those that are Rayleigh-stable by a 
sufficiently small margin can become turbulent through a nonlinear instability associ¬ 
ated with a subcritical bifurcation from infinite Reynolds number. In the non-rotating 
case the laminar state admits algebraically growing infinitesimal disturbances that are 
damped only on the viscous timescale Re/S. The nonlinear terms of the model allow 
perturbations of finite amplitude to be sustained and the system evolves to a non-trivial 
state of statistically steady turbulence. This behaviour of the model strongly resembles 
the theory of subcritical transition to turbulence, developed by Trefethen et al. (1993) 
and others, involving the transient amplification of disturbances by a non-normal oper¬ 
ator and a cooperative nonlinear feedback. The closure model that we work with has 
the advantage of being able to represent the final turbulent outcome of the transition 
process. 

The analysis of Reynolds-stress models in homogeneous shear flow in terms of a non¬ 
linear dynamical system is not unique to our work (see, e.g., Speziale, Gatski & Mac 
Giolla Mhuiris 1990). However, the simplicity of our model permits an exhaustive study 
of its dynamical properties and, by including the effects of a finite Reynolds number, 
we are able to make a connection with the theory of subcritical transition to turbulence 
in which the laminar state has a basin of attraction that diminishes as the Reynolds 
number is increased. Similar techniques of analysis could of course be applied to more 
sophisticated closure models and we believe that our findings are to some extent generic. 

The turbulent solutions are anisotropic as a result of shear and rotation, and in the limit 
of large Reynolds number the shear stress behaves as in Prandtl’s mixing-length theory, 
but with a prefactor that depends on the Rayleigh discriminant (see equation (13.1511 1. 
As such, our model naturally captures the reduction, and eventually suppression, of the 
turbulent energy dissipation for rapidly rotating flows (Speziale et al. 1998). 

When applied to wall-bounded turbulent shear flows fSection rOll . the model predicts 
the occurrence of a universal velocity profile close to a wall at large Reynolds number. 
Outside the viscous sublayer, this profile has the logarithmic form predicted by Prandtl’s 
mixing-length theory, and we derive two accurate constraints on the three parameters 
from matching the most recent experimental data (Zagarola & Smits 1998). 

We have also investigated in some detail the predictions of the model for the occurrence 
and the characteristics of turbulent states in Couette-Taylor flow without end-effects 
(Section 0J. Here, depending on the ratios of the radii and angular velocities of the two 
cylinders, the distribution of the Rayleigh discriminant of the laminar solution may be 
such that a local analysis would predict either linear instability or nonlinear instability 
or complete stability in different regions of the flow. Furthermore, once turbulence sets 
in, the angular velocity distribution and the corresponding Rayleigh discriminant are 
significantly modified from those of the laminar solution. Therefore a wide variety of 
behaviour is possible, including the existence of mixed states, in which the turbulence 
is localized. It is worth noting that although we have restricted the present analysis to 
solutions of maximal symmetry in the Couette-Taylor system, our model may admit 
classes of more general solutions. For instance, relaxing the assumption of azimuthal and 
axial translational symmetry could in principle help explain observed phenomena such 
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as spiral turbulence (Coles 1965; Hegseth et al. 1989) in which regions of laminar and 
turbulent flows coexist separated by a helical interface. 

By fitting the remaining parameter of our model, we are able to account quite well 
for the qualitative behaviour and quantitative torque measurements in historical exper¬ 
iments on Couette-Taylor flow by Wendt (1933) and Taylor (1936a), which have not 
been superseded. As an unexpected bonus, the model captures reasonably accurately the 
appearance of Taylor vortices at the onset of linear instability. 

It is appropriate to discuss at this point some of the implications of the parameteriza¬ 
tion used in our model. The ratio C' 2 /C'i represents the propensity of the turbulence to 
return to isotropy. In a local analysis, it is found to determine the critical value of the 
dimensionless Rayleigh discriminant for which high-Reynolds-number turbulence can be 
sustained. Why these properties should be related can be explained with reference to the 
system of equations (TOt . When the Rayleigh discriminant = 211(217 — S) is positive, 
the terms AflR^y and 2(5' — 2Vt)Rxy in the equations for Rxx and Ryy have opposite 
signs. Therefore either Rxx or Ryy lacks a positive source, and the turbulence must de¬ 
cay, unless the C 2 term comes into play. For a given positive Rayleigh discriminant, the 
isotropizing tendency must be sufliciently great if the turbulence is to be sustained. 

In this work, although we have proposed values of the coefficients Ci , C 2 and after 
fitting experimental data, these values are only tentative and approximate and we do not 
claim that such a simple model can provide great quantitative accuracy in comparison 
with currently available closure models (cf. Choi & Lumley 2001) . In particular, the 
ratio C' 2 /C'i is only weakly constrained through a comparison with Wendt’s (1933) data, 
with a wide plausible range of roughly 0.4 to 1. 

We draw attention again to the important problem of the hydrodynamic stability of 
circular Keplerian motion in astrophysical accretion discs, in which the angular velocity 
profile 17 oc r~^l^ is enforced by gravitational dynamics, not through the boundaries. 
While Richard & Zahn (1999) sought to apply the findings of Couette-Taylor experi¬ 
ments to accretion discs, our investigation of turbulent Couette-Taylor flows suggests 
that caution is required in making such associations. According to our model (taking 
C' 2 /C'i < 8/3), Keplerian rotation is likely not to support statistically steady turbulence 
in a local analysis, and may be nonlinearly stable no matter how large the Reynolds num¬ 
ber. We also find that this does not contradict the experimental finding that Couette- 
Taylor flow with a stationary inner cylinder becomes turbulent at large Reynolds number, 
and is even consistent with the possibility that a wide-gap Couette-Taylor flow with the 
cylinders in a Keplerian ratio may be turbulent. In addition, Couette-Taylor experiments 
are always contaminated to some degree by end-effects. We suggest that Couette-Taylor 
experiments may be of limited applicability to the study of the nonlinear stability of 
Keplerian rotation, and that it can be instead most usefully addressed in local numerical 
models such as that of the shearing box (e.g. Balbus & Hawley 1998), which are free from 
end-effects and also from the type of radial boundary conditions that induce boundary 
layers. To date, no instability has been found in such models, and it would be valuable 
to test this finding to very high Reynolds numbers. 

Finally, we emphasize that the philosophy behind the construction of simple turbulence 
models such as the one adopted here is applicable to a range of more complex problems 
such as magnetohydrodynamic turbulence (Ogilvie 2003), convection, or mixing in strat¬ 
ified shear flows. Such extension is the subject of current investigations. 
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of the Royal Society through a University Research Fellowship. 
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Appendix A. Model equations in a general orthogonal curvilinear 
coordinate system 

Using Batchelor’s (1967) notation, the equivalent of equation 12.811 for the evolution of 
the Reynolds stress tensor in general orthogonal curvilinear coordinates is 


dt 
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+ - 


hj^ dxiz 
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■ij ^ dui 
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(Al) 


where {hi, /i 2 , / 13 ) is (1,1,1) for Cartesian coordinates {x, y, z), or (1, r, 1) for cylindrical 
coordinates {r,(p,z). 

In the Cartesian case, the viscous correction terms follow from the decomposition 


viu'i^kku'j + m'- fcfcw') = v{{u{uj),kk - 2?x'_fcu'-fe) , 


(A 2) 


where the first term on the right-hand side describes a viscous diffusion of the stresses 
and the second term describes a direct decay of the stresses, which we then model as 

-a^i?y/L 2 . 

In order to obtain the form for the viscous corrections in general curvilinear coordi¬ 
nates, we follow the same method used in the Cartesian case: we isolate from the original 
terms v(y^u')iu'^ + v{V'^u')ju{ the tensor decay term —2v{yu'V"^u')ij (which is the 
covariant equivalent of —2vu{ k'^j k)^ where Vm' is the matrix defined by its columns 
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In this expression, u' = u^ei -|- 1^262 + 1(363 and the derivatives of the unit vectors 
(ei, 02 , 63 ) are given by Batchelor (1967, p. 598). We model the decay terms as —CyvRij/L'^, 
and keep the remaining diffusion terms unchanged thereby defining the Laplacian of the 
second rank tensor (V^i?)ij. Hence with this method the viscous terms in equation 12.911 
are then simply 

-^%+KV"R)y, (A 4) 
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in cylindrical coordinates for instance. 

For an axisymmetric flow with translational symmetry in z and u = rf2(r) in a 
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cylindrical geometry the evolution equations for the Reynolds stress tensor are 




dt 
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Appendix B. Asymptotic solution of the turbulent Couette—Taylor 
flow for large Reynolds number and large Rossby 
number 

In the limits of large Rossby number (as expected for a small gap) and large Reynolds 
number, the angular momentum equation can be approximated (to first order in 
Ro~^) by 

r^K^L^S\S\ + 2r^K^L^n\S \^ , (B 1) 

C^2 2i7Thp 

where k is the von Karman constant defined in equation (Tna . This provides a quadratic 
equation for S which can be inverted, and yields (in the same limit) 


S = 
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This equation must be integrated separately in each intervals (ri,rm] and [r,„,ro): in 
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Let a = 9 C 1 IC 2 + 2. The integrals in equations fm and (IB 4ll have a logarithmic 
singularity as r ^ and r ^ ro, which can be isolated as 



dr' = f{r; r^, r*, a) - “ In 



(B5) 


This defines the function / uniquely. The logarithmic singularity naturally matches onto 
the boundary layer solutions near the walls. 
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We use the results of Section E21 to write the boundary layer solution explicitly. Near 


r = Ti, but outside the viscous sublayer, 




2TThp 




(B6) 


whereas near r = r„ 


O. ^ O , sign(5') / \T\ 
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Matching the inner (IFHll . ETIl and outer iTBlll . (TBUl solutions near the walls, and 
continuity across the mid-point provides an equation for the torque, which extends 
the Prandtl-von Karman skin-friction law for Couette-Taylor flows: 
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where we have definedf Re = — floK^’o ~ C/ = T/Re^phv'^ and the friction 

law coefficients M and N as 
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Note that in the limit where the contribution from rotation on the Reynolds stresses is 
neglected {Ro~^ = 0), the solution in the bulk of the fluid can be written out as equations 
(IB 3ll and IjB 4ll with a = 2 and 


-I-In r —In Cm J ■ (BIO) 

In Figure CHI we compare the velocity profiles obtained by numerical integration to 
those derived from asymptotic analysis, for three different gap widths. We find that the 
asymptotics only provide accurate results for d/vo < 0.02. This somewhat disappointing 
range of applicability of the Ro ^ 1 asymptotic analysis is due to the great efficiency 
with which turbulence redistributes the shear, which reduces the Rossby number in the 
interior of the flow compared to that of the laminar solution. 


fir;rm,ri,2) = 
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